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We present Monte Carlo simulations for a classical antiferromagnetic Heisenberg model with both 
nearest (Ji) and next-nearest (J2) exchange couplings on the square lattice in the presence of non- 
magnetic impurities. We show that the order-by-disorder entropy selection, associated with the 
Ising-like phase transition that appears for J2/J1 > 1/2 in the pure spin model, is quenched at low 
temperature due to the presence of non-magnetic impurities. Evidences that a new competing order 
is stabilized around the impurities, and in turn induces a re-entrance phase transition are reported. 
Implications for local magnetic measurement of the parent compound of iron pnictides are briefly 
discussed. 



Unconventional superconductivity occurs in the prox- 
imity of magnetically ordered state in many materials 
PQ [2]. Understanding the magnetic phase of the par- 
ent compound is an important step towards understand- 
ing the mechanism of superconductivity. Unlike the case 
of the cuprates, magnetism and its underlying electronic 
state in the iron pnictide superconductor BaFe2As2 [3] is 
still not well understood. Many low-energy probes such 
as transport [4 j, scanning tunnelling microscopy [5] and 
angle-resolved photoemission spectroscopy [6] have mea- 
sured strong in-plane anisotropy of the electronic states, 
but there is no consensus on its physical origin. 

It was suggested from first principle calculations [7 
that the origin stems from orbital order, although the 
obtained anisotropy in the resistivity is opposite to the 
one found experimentally [8 . A more likely scenario is 
related to spin density wave instabilities, which is sup- 
ported by recent neutron diffraction measurements [9], 
and stems from Fermi surface nesting of electron and 
hole pockets. In the latter picture, the nematic magnetic 
order introduce in turn an orbital polarization, since the 
electron pockets at Q = (7r, 0) and Q = (0, tt) have d yz 
and d xz orbital characters respectively p~0j[TT]. The ex- 
act nature of the magnetic ground state remains however 
unclear. It was suggested that the commensurate AFM 
order can also be described within a local moment picture 
that may become relevant in the presence of moderately 
large electronic correlations and can be quantified, for 
example, in terms of the Heisenberg model with both 
nearest- (Ji) and next-nearest (J2) exchange couplings 
in the frustrated regime J 2 > 2Ji [HI- The first esti- 
mation of the coupling constants Ji 5 2, obtained by fit- 
ting the experimental spin density wave excitation spec- 
tra, yielded parameters which are not in the frustrated 
regime [13]. However, it was later shown that the fits 
of the experimental data included energy scales beyond 
lOOmeV, which are not well described by magnon excita- 
tions [14 . A more careful study, including the itinerant 
character of the electrons, suggested that indeed the pnic- 
tides are in the frustrated regime [15] . This scenario was 
also recently also supported by first-principle calculations 



for selenium based compounds (KFe2Se2 ) [16 , and the 
iron pnictides BaFe2As2 and KFe2Se2 compounds were 
proposed as experimental realisations of a layered J1—J2 
spin model in the collinear regime [17 . In particular, it 
has been suggested both experimentally [18] and theoret- 
ically [19] that non-magnetic impurities have a dramatic 
impact on the magnetic and superconducting properties. 

In this Letter, we address the question of the interplay 
between the frustration, induced by the exchange cou- 
pling, and the disorder induced by the imperfections of 
the crystallographic structure. Firstly, we consider non- 
magnetic impurities, and then extend the calculations 
to magnetic impurities. Since density functional calcu- 
lations, and quite generally quantum based calculations, 
are limited to relatively small unit-cell and cannot tackle 
the issue of large super-cell structures, we limit our cal- 
culations to the classical case, and carry out Monte Carlo 
calculations of the J1 — J2 model in the presence of non- 
magnetic impurities. The methodology and implementa- 
tion were discussed in Ref. [2Q1 [21] , a short summary is 
given hereafter. 

The Heisenberg hamiltonian reads: 

H = y £ t J 1 S i >S j + J 2S,-S„ (1) 

(iJ) <W» 

where are 0(3) spins on a periodic square lattice with 
N = L x L sites, (i, j) and indicate the sum 

over nearest and next-nearest neighbors, respectively. J\ 
sets the energy scale, and J2/J1 = 0.55 is used through- 
out the rest of the paper. In absence of disorder, the 
magnetic vector is Q = (7r, tt) for J2/J1 < 0.5, and for 
J2/J1 > 0.5 the ground state is continuously degener- 
ate, but the entropy selection reduces the 0(3) symmetry 
of the ground state to Z2 at finite temperature, select- 
ing the states with antiferromagnetic spin correlations in 
one spatial direction and ferromagnetic correlations in 
the other (Q = (0,7r), (7r,0)). This is the so-called order 
by disorder scenario, and the associated discrete symme- 
try breaking drives a finite temperature Ising-like phase 
transition [17]. To characterize this transition, it is use- 
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FIG. 1: (Color online) a) Typical spin configuration obtained 
at T — 10 _6 Ji. The filled circle indicates the location of the 
impurity, and the rectangle highlights the region around the 
impurity where the spins deviate from the 90° ordered state 
for J2/J1 = 0.55. b) Angle between the spins connected to 
the impurity with the horizontal axis (angle a) as obtained 
from the Monte Carlo (squares) and compared with a sim- 
ple variational criteria where only the spins connected to the 
impurity are tilted (line). Inset: distortion angle a as a func- 
tion of the distance from the impurity (in units of the lattice 
spacing ao) obtained by Monte Carlo for J2/J1 = 0.55. 



FIG. 2: (Color online) Spatially resolved order parameter 
a) for the 90° spin order (M 90 (x)) and b) for the Ising order 
(M2(x)) obtained at T = 10 _3 Ji. The respective order pa- 
rameters evaluated at a lower temperature T = 10 _4 Ji are 
shown in c) and d). In all the calculations above the cluster 
contain LxL=120xl20 sites and the impurity is located 
at the center of the cluster at (x,y) = (60, 60). e) Color plot 
of Mqo(x) at T — 10 _3 Ji (blue and red are respectively the 
minimum and maximum). All calculations above are done for 
J 2 /J x = 0.55. 



fill to construct, from the original spin variables S^, an 
effective Ising variable on the dual lattice: 

M 2 {x) = (S i -S k )-(S j -S l ), (2) 

where (i,j,k,l) are the corners with diagonal (i,k) and 
(j, I) of the plaquette centered at the site x of the dual lat- 
tice, and we define its normalized counterpart as Z 2 {x) = 
M2 (x) 1 1 M2 (x) I . In this way, the two collinear states with 
Q = (7r, 0) and Q = (0,7r) can be distinguished by the 
value of the Ising variable, Z 2 (x) = ±1. 

We first discuss the symmetry of the magnetic order 
for the case of a single impurity. The ground state for 
J2/J1 > 0.5 in the absence of impurity is continuously 
degenerate and is characterised by a bi-partite lattice, 
with two distinct anti-ferromagnetically ordered states 
on respectively each sub-lattice, and is the angle be- 
tween the two magnetic directions. Our Monte Carlo 
calculations show that introducing a single impurity lifts 
the former continuous degeneracy, and selects the state 
with 6 = 90°, as shown in Fig. [l]a, and in agreement 
with a prediction by Chris Henley [22], who suggested the 
name anticollinear to describe the state with 6 = 90°. 
This state was also suggested as a stable phase of ferro- 
pnictides [T5] . 

This selection is induced by a local energy optimiza- 
tion around the impurity site. Indeed, the spins on the 
impurity's nearest neighbor sites slightly distort in order 
to align ferromagnetically, in order to optimize locally 
the energy with their own neighbors once an impurity is 
introduced into the 90° state. We show in Fig. [l]b the 
obtained deviation of the spins connected to the impurity 
from the bulk 90° magnetic state. We obtain, as shown in 
the inset of Fig.[l]b, that the bonds affected by this local 



distortion are in the very near vicinity of the impurity. 
Noteworthy, treating the distortion angle a as a simple 
variational parameter, and neglecting the distortion of 
the spins which are not nearest neighbors of the impurity, 
lead to a very good estimate of the Monte Carlo result. 
The good agreement between the two methods confirms 
that the energy optimization around a single impurity is 
essentially local in space, and thus the selection of the 90° 
magnetic state is driven by a local energetic optimization 
process, in contrast to the order by disorder selection of 
the Q = (0, 7r), (7r, 0) states, which is an entropic selec- 
tion. 

Next, we investigate the relevance of this state as 
a function of temperature and impurity concentra- 
tion. This anticollinear order is characterized similarly: 
Mqq(x) = I (Si — Sfc) A(Sj — Si)\. Interestingly, the devia- 
tion is large for J 2 ~ ^i/2, and decreases significantly as 
J 2 increases. The latter is easily explained since, when J 2 
become large, the four nearest neighbour spin of the im- 
purity get antiferromagnatically aligned, and hence the 
distortion a is not energetically favored. Interestingly, 
the energy optimization AE = E(a) — E(a = 0) ~ 
— 0.7Ji is quite significant around J2/J1 = 0.5 and hence 
we expect the 90° magnetic phase to be a relevant mag- 
netic phase for pnictides, since first-principle calculations 
obtain values for J2/J1 close to 0.5 (such as J2/J1 ~ 0.75 
for KFe2Se2 PS]). We now turn to the discussion of the 
single impurity problem at small but finite temperature. 
The 90° spin order around a single impurity is shown 
in Fig. [2]a,c. Indeed, the 90° order is not breaking a 
discrete symmetry, such as the Ising symmetry broken 
by the Z2 Ising parameter defined in equation ([2|, and 
is not a stable thermodynamic phase. Thus, we observe 
that the 90° order does not develop long-range correla- 
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tions but is rather stabilized around the impurity within 
a finite region, which defines its correlation distance. At 
finite temperature, there is a competition in the free en- 
ergy F = E — TS between, on the one hand, the local 
energy optimization in the vicinity of the impurity which 
favors the 90° spin order, and on the other hand the en- 
tropy selection which favors the Ising Q = (0, 7r), (tt, 0) 
states. Notwithstanding that the 90° state is energeti- 
cally favored and is stabilized at short distances from the 
impurity, we find that thermal fluctuations screen the 
impurity at long distances such that the system recovers 
the entropically selected Ising states far from the impu- 
rity. Hence, there are mainly three different spatial scales 
in this problem: i) The direct proximity of the impurity, 
where the 90° state distorts to optimize the energy (£i), 
ii) the correlation length of the 90° state for which this or- 
der survives the entropic selection (£2 ) , hi) and finally the 
correlation length of the Ising state (£3). Furthermore, 
the interplay between the spatial scales is temperature 
dependent, in particular the size of the 90° cluster, as 
seen from Fig. |2ja obtained at T = 10 -3 J\ and Fig. [2]c 
obtained at a smaller temperature T = 10 _4 Ji. In addi- 
tion, we note that the shape of the 90° cluster is highly 
asymmetric. To understand this point, we measured the 
Ising order parameter M2 in the same calculation (see 
Fig. |2jb,d), and as expected we observe a concomitant 
reduction of the Ising order in the region where the 90° 
order is large. Interestingly, we find that the shape of 
the 90° cluster correlates with the Ising order: the lat- 
ter is an ellipsoid (the large axis is denoted by Ti and 
the small axis by I^), we find that Ti is along the e y 
(e x ) direction for the corresponding Z2 = +1 (Z2 = — 1) 
Ising state. In particular, as shown in Fig. [2]e, the spins 
along Ti (r 2 ), highlighted by the white stripes, corre- 
sponds to the ferromagnetically (anti-ferromagnetically ) 
aligned spins of the colinear Ising states. Hence, due to 
the antiferromagnetic J\ coupling, the deviations from 
the pure Ising state are energetically favorable along Y\ , 
and costly along T2, which explains why the screening of 
the impurity is more effective in one direction than the 
other. 

Let us now extend the discussion to a finite concentra- 
tion of impurities 5 (see Fig. [3ja). We carried out Monte 
Carlo calculations for a L x L = 40 x 40 cluster with dif- 
ferent magnetic impurities dilutions and we averaged 
the physical observables over 32 random configurations of 
impurities. We find for small concentration 5 = 0.125% 
only a weak effect on the Ising order. In particular, we 
observe the Ising-like cross-over at T w 0.2Ji, as shown 
by the sharp drop of the order parameter as this temper- 
ature, and a very small dip in the Z 2 order parameter 
at T « 0.005Ji. We note that transitions belonging to 
the Ising universality class {y = 1 and the dimension 
d = 2) do not satisfy the well known Harris criteria [23] , 
which assess that phase transition with vd > 2 are un- 
affected by the disorder. However, as numerical studies 




FIG. 3: (Color online) a) Temperature dependence of the 
spatially averaged Ising order parameter Z2 for various impu- 
rity concentrations 5 obtained for a L x L = 40 x 40 lattice. 
The Ising order is suppressed at low temperature by the pres- 
ence of impurities, and at large temperature by strong ther- 
mal fluctuations. The dashed lines are guide to the eyes to 
track the Ising crossovers at low and high temperatures, b) 
Corresponding phase diagram in impurity density 6 and tem- 
perature T. All calculations above are done for J2/J1 = 0.55 
and were averaged over 32 random configurations of disorder. 



show that the 2D Ising model is weakly affected by the 
disorder [24 , our results show that the cross-over tem- 
perature at T « 0.2 J 1 is not strongly affected by moder- 
ate disorder, and hence suggest that the phase transition 
survives in the presence of impurities. A formal proof 
would require a detailed finite size scaling analysis and 
goes beyond the scope of this work. 

At larger concentrations 0.125% < 5 < 11%, we find 
a steady decrease of the Ising order at small tempera- 
tures, as highlighted by the dashed lines in Fig. [3Ja. For 
instance, at S = 2% we observe a quench of the Ising 
order for T < 0.025Ji, and the entropic selection kicks 
in for 0.025 < T < 0.19 where we obtain the Ising or- 
dered phase, and finally at higher temperatures t > 0.19 
the system is a paramagnet. This scenario is commonly 
called re- entrance phase transition. We note that the 
impurities mainly affect the Ising order at low temper- 
atures, and the Ising-like transition near T « 0.2Ji is 
moderately affected by impurities at small and mod- 
erate dilutions. Beyond a critical dilution 5 C « 20%, 
we do not observe the presence of the collinear or anti- 
collinear states. Indeed, in two dimensions and for an 
impurity dilution S c = 1/9, there is one impurity in av- 
erage connected to every spin of the lattice, so that the 
local distortions and subsequent local energy optimiza- 
tions start to prevail over both the Ising phase and the 
90° local spin order. The phase diagram is summarized 
in Fig.[3]b. Note that the reentrant behavior of the Ising 
phase agrees with the prediction of Ref. [22]. However, 
the rest of the phase diagram of Ref. [22] cannot be com- 
pared to the present results. Indeed, in the case of the 
XY model studied in Ref. [22 , the local chirality of the 
anticollinear order defines an Ising variable, and the anti- 
collinear phase must be separated from the paramagnetic 
phase by a phase transition. By contrast, no long-range 
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FIG. 4: (Color online) Real space map of M 90 obtained at 
a) T = 10" 3 Ji, b) T = 10" 4 Ji and c) T = 10" 5 Ji for a given 
configuration of 16 non-magnetic impurities in a L x L — 
120 x 120 lattice (8 = 0.1%). The colours range from black 
(minimum) to red (maximum), d) Critical dilution S c as a 
function of the temperature obtained by the Monte Carlo data 
of Fig.[3]a (triangles), and by a simple criteria comparing the 
length of the short (dark gray circles) and long (light gray 
circles) ellipsoidal axis of the 90° order cluster (see discussion 
in the text), e) Ti (the long axis of the 90° order ellipsoid) 
in units of ao as a function of r = Si mp /Siat, where Simp is 
the spin of the magnetic impurity, and Si a t is the spin of the 
correlated element of the compound. All calculations above 
are done for J2/J1 = 0.55. 



order can exist for the vector chirality at finite tempera- 
ture for Heisenberg spins, and the low-temperature phase 
below the Ising phase can be smoothly connected to the 
paramagnetic phase. To study this crossover would re- 
quire to study larger dilutions and to perform a system- 
atic finite size scalings, which goes beyond the scope of 
the present work. 

Interestingly enough, the re-entrance phenomena ob- 
served in Fig.[3]a can be explained at small impurity con- 
centration on the basis of the single impurity results at 
finite temperature. Let us consider a given configuration 
of impurities for a finite impurity dilution of S = 0.01 at 
T = 10- 3 , 10- 4 , 10- 5 Ji (see Fig.[Z}a,b,c). The 90° order 
stays localized around the impurities at high temperature 
(Fig. [4ja), but forms superstructures connecting the im- 
purities (Fig. |4]b) upon lowering the temperature until 
it finally spreads through the whole lattice (Fig. |4jc). 
This process is very similar to a percolation transition, 
and can be captured within a very simple argument: The 
90° order become long-range when the 90° order corre- 
lation length £90 (T) is of the order of the mean distance 
A between the impurities. We note that the 90° is an 
ellipsoid, and hence we obtain a lower and upper bound 
on the critical dilution, by comparing the short and long 
axis respectively to the mean impurity-impurity distance. 
Indeed, we compare in Fig. |4jd the critical dilution ob- 
tained from the Monte Carlo data of Fig.[3ja and with the 
critical dilution obtained by this simple argument. We 
find that this argument provides a reliable estimate of 
the critical dilution. This confirms that the re-entrance 



of the Ising order is driven by the competition of differ- 
ent length scales associated to the Ising and 90° order 
parameters. 

Finally, we generalized our calculations to magnetic 
impurities with non zero spin (see Fig. [4je). Remark- 
ably, we find that the 90° order cluster around a sin- 
gle impurity is not significantly affected by the spin of 
the impurity Si mp , as long as the ratio of the spin of 
the impurity to the magnetic element of the compound 
r — Simp/ Si a t remains smaller than w 0.6. This suggests 
that the re-entrance phase transition does not strictly re- 
quire non-magnetic impurities but could also be present 
for instance in the case of Ni impurities in BaFe2As2 , 
where the spin of Ni is « 40% of the spin of Fe[25j. So 
it will be very interesting to see if the order proposed in 
this paper can lead to an alternative interpretation of the 
NMR results in BaFe2As2 and maybe help clarifying the 
origin of the lineshapes in Ni and Zn doped samples [25 . 

In conclusion, we have carried out a systematic study 
of the effect of non-magnetic impurities in a frustrated 
Heisenberg model. We reported that for J2/J1 > 0.5 the 
continuous degeneracy of the ground state induced by 
the frustration is lifted due to a local optimization in the 
vicinity of a single impurity. The energy gain favors an- 
ticollinear order, which consists of bipartite lattices sup- 
porting Neel states entangled with a 90° angle. This or- 
der, energetically favored, competes with the Ising order, 
entropically favored, and at long distance we find that 
the impurity is screened by thermal fluctuations. This 
results in a rich phase diagram with the stabilization of 
anticollinear order at low temperature as soon as a finite 
concentration of impurities is present followed by a reen- 
trant collinear phase upon increasing the temperature. 
Moreover, we have shown that the structure around the 
impurity locally departs from the purely anticollinear or- 
der. This effect is large when J2/J1 is close to 1/2, as in 
the pnictides, and should be detectable by local probes 
such as NMR. 
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